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Abstract 

Self-calibration  techniques  have  been  used  extensively  in  co-ordinate  metrology.  At  their 
most  developed,  they  are  able  to  extract  all  systematic  error  behaviour  associated  with 
the  measuring  instrument  as  well  as  determining  the  geometry  of  the  artefact  being 
measured.  However,  this  is  generally  at  the  expense  of  introducing  extra  parameters 
leading  to  moderately  large  observation  matrices.  Fortunately,  these  matrices  tend  to 
have  sparse,  block  structure  in  which  the  nonzero  elements  are  confined  to  much  smaller 
submatrices.  This  structure  can  be  exploited  either  in  direct  approaches  in  which  QR 
factorisations  are  performed  or  in  iterative  algorithms  which  depend  on  matrix-vector 
multiplications.  In  this  paper,  we  describe  self-calibration  approaches  associated  with  high 
accuracy,  dimensional  assessment  by  co-ordinate  measuring  systems,  highlighting  how  the 
associated  optimisation  problems  can  be  presented  compactly  and  solved  efficiently.  The 
self-calibration  techniques  lead  to  uncertainties  significantly  smaller  than  can  be  expected 
from  standard  methods. 

1 Introduction 

An  important  activity  in  metrology  is  the  calibration  of  instruments  and  artefacts.  Cal- 
ibration defines  a rule  which  converts  the  values  output  by  the  instrument’s  sensor(s) 
to  values  that  can  be  related  to  the  appropriate  standard  (SI  or  derived)  units.  Import- 
antly, to  these  calibrated  values  it  is  required  to  assign  uncertainties  that  reliably  take 
into  account  the  uncertainties  of  all  quantities  that  have  an  influence.  As  a consequence, 
the  size  and  complexity  of  the  computational  tasks  associated  with  the  data  analysis  can 
be  significant,  even  for  instruments  that  appear  to  be  of  simple  design  and  operation. 
It  is  thus  beneficial  to  design  and  implement,  algorithms  that  are  efficient,  with  respect 
to  computation  and  memory.  Fortunately,  many  of  the  calibration  problems  give  rise  to 
systems  of  equations  with  a well  defined  sparsity  structure. 

The  rest  of  this  paper  is  organised  as  follows.  In  Section  2 we  review  least  squares 
approaches  to  calibration  problems  and  go  on  to  describe  self-calibration  problems  in 
co-ordinate  metrology  in  Section  3.  Sections  4 and  5 describe  solution  methods  for  two 
types  of  sparsity  structure.  Our  concluding  remarks  are  given  in  Section  6. 

2 Least  squares  solution  to  calibration  problems 

In  many  calibration  problems,  the  observation  equations  involving  measurements  y. 
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can  be  expressed  as  y*  = <pi  (a)  + e,,  where  4>t  is  a function  depending  on  parameters 
a = (oi, . . . , a.„)T  specifying  the  behaviour  of  the  instrument,  and  e,  represents  random 
measurement  error.  For  a set  of  measurement  data  {y,} j",  best  estimates  a*  of  the 
calibration  parameters  a are  determined  by  solving 

m 

min  ^2ff( (2.1) 
i=i 

where  /*( a)  = y»  — <^>,(a).  The  most  common  approach  to  solving  this  problem  is  derived 
from  the  Gauss-Newton  algorithm;  see,  for  example,  [5].  If  a is  an  estimate  of  the  solution 
and  J is  the  Jacobian  matrix  defined  at  a by  Jij  = dfi/daj , then  an  updated  estimate 
of  the  solution  is  a + p,  where  p solves  the  Jacobian  system 

Jp  = -f, 

in  the  least  squares  sense.  Starting  with  an  appropriate  initial  estimate  of  a,  these  steps 
are  repeated  until  convergence  criteria  are  met. 

A numerically  stable  method  of  solving  the  Jacobian  system  is  to  find  a factorisation 
J = QR,  where  Q is  an  m x n orthogonal  matrix  and  R is  an  upper-triangular  matrix 
of  order  n (see,  e.g.,  [1,  6]).  The  solution  p is  determined  efficiently  by  solving  the 
upper-triangular  system 

Rp  = -Qt  f, 

using  back  substitution.  The  matrix  Q can  be  constructed  using  either  Householder 
reflections,  which  process  the  Jacobian  matrix  a column  at  a time,  or  Givens  plane 
rotations,  which  process  the  matrix  row-wise.  For  either  approach  the  orthogonal  fac- 
torisation requires  0(mn2)  operations. 

An  alternative  to  the  direct  approaches  to  solve  matrix  equations  is  to  use  iterative 
procedures  based  on  conjugate  gradients.  The  advantage  of  these  approaches  is  that  they 
involve  only  matrix-vector  multiplications  and  for  sparse  matrices  these  multiplications 
can  be  made  efficient.  In  particular,  the  LSQR  algorithm  of  Paige  and  Saunders  [7] 
implements  an  iterative  approach  to  solving  linear  least  squares  problems. 

Often,  linear  equality  constraints  on  the  parameters  of  the  form  CTa  = c,  where  C 
is  an  n x p matrix,  p < n,  are  required  to  eliminate  degrees  of  freedom  in  the  problem. 
However,  we  can  use  orthogonal  projections  to  eliminate  these  constraints.  Suppose  C 
is  of  full  column  rank  and  has  QR  factorisation 

C=[Vl  V2] 

where  Vj  and  V2,  respectively,  are  the  first  p and  last  n — p columns  of  the  orthogonal 
factor  V.  If  ao  is  a solution  of  CTa  = c,  then  for  any  (n  — p) -vector  a,  a = a0  + Vffi 
automatically  satisfies  the  constraints  and  the  optimisation  problem  can  be  reformulated 
as  the  unconstrained  non-linear  least  squares  problem 

m 

min  Vfyfyao  + V2a), 

a z ' 

7=1 


S 
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involving  the  reduced  set  of  parameters  a.  We  note  that  the  associated  Jacobian  matrix 
is  simply  J = JV 2,  where  J,j  = dfi/daj,  as  before. 

Unfortunately,  even  if  J has  structure  J = JV2  could  be  full.  For  indirect  approaches, 
this  is  of  little  consequence  since  the  matrix-vector  multiplications  can  be  formed  in  two 
stages  (e.g.,  y = V2x,  z = Jy)  each  of  which  can  be  implemented  efficiently.  For  a direct 
approach,  it  may  be  possible  to  implement  the  constraints  in  such  a way  as  to  minimise 
the  amount  of  fill-in  during  the  orthogonal  factorisation  stage. 

3 Self-calibration  problems  in  co-ordinate  metrology 

Co-ordinate  metrology  is  concerned  with  defining  the  geometry  of  two  and  three  dimen- 
sional artefacts  from  measurements  of  the  co-ordinates  of  points  related  to  the  surface 
of  the  artefacts.  It  is  a key  discipline  in  quality  and  process  control  in  manufacturing 
industy.  In  a (conventional)  co-ordinate  measuring  machine  (CMM)  with  three  mutu- 
ally orthogonal  linear  axes,  the  position  of  the  probe  tip  centre  is  inferred  from  scale 
readings  on  each  of  the  three  machine  axes.  In  practice,  CMMs  have  imperfect  geometry 
with  respect  to  the  straightness  of  the  axes,  the  squareness  of  pairs  of  axes  and  rotations 
describing  roll,  pitch  and  yaw,  and  these  systematic  errors  have  to  be  taken  into  account 
if  the  accuracy  potential  of  the  CMM  is  to  be  more  fully  realised.  Two  approaches  can 
be  adopted  to  nullify  the  effect  of  these  systematic  errors.  The  first  - error  mapping  - 
involves  performing  a set  of  experiments  to  characterise  as  completely  as  possible  the 
error  behaviour  of  the  instrument  and  then  use  error  correction  software  to  produce  more 
accurate  co-ordinate  estimates.  The  disadvantages  of  this  approach  are,  firstly,  the  set 
of  experiments  is  expensive  to  perform  and,  secondly  and  more  importantly,  the  error 
behaviour  of  the  CMM  is  likely  to  drift  so  that,  for  example,  an  error  correction  valid  on 
Monday  will  only  be  partially  valid  on  Friday  and  may  be  of  limited  value  a month  later. 
The  second  approach  - self-calibration  - attempts  to  use  any  approximate  symmetries, 
rotational  or  translational,  of  the  artefact  so  that  systematic  errors  associated  with  the 
measuring  system  are  identified  as  part  of  the  measurement  process  [4] . The  advantage 
of  this  method  is  that  the  effect  of  systematic  error  behaviour  of  the  instrument  is  can- 
celled out  arid  the  accuracy  of  the  measurements  are  limited  only  by  the  smaller,  random 
component. 

3.1  Calibration  of  reference  artefacts  in  2-dimensions 

As  an  example,  we  consider  the  accurate  calibration  of  2-dimensional  artefacts  by  a two 
dimensional  CMM.  The  artefacts  define  the  location  of  targets  nominally  aligned  in  a 
grid  pattern.  Let  y j,  j = 1, . . . ,ny,  be  the  locations  of  the  targets  in  a fixed  frame  of 
reference,  and  let 

yj,k=T{yj,tk) 

be  the  location  of  the  jth  target  in  the  fcth  measuring  position.  Here,  the  roto-translation 
T is  specified  by  three  parameters  t defining  the  translation  vector  and  angle  of  rotation. 

We  suppose  the  systematic  error  of  the  two  dimensional  CMM  can  be  expressed  as 

x*  = x*(x,  b)  = x 4-  e(x,  b), 
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where  x*  are  the  true  point  co-ordinates,  x are  the  indicated  point  co-ordinates  output  by 
the  machine  and  e(x,  b)  is  the  error  correction  term  depending  on  x and  error  parameters 
b.  For  instance,  suppose  the  model  describes  scale  and  orthogonality  errors  so  that 

x*  = x(l  + 61)  + 2/(1  + 62)  sin 63,  y*  = y(l  + b2)  cos 63. 

If  Xj  is  the  measurement  of  the  jth  target  with  the  artefact  in  the  ktb.  position  then  the 
associated  observation  equation  is 

Xi  + e(xj,b)  = yjtk  + e{.  (3.1) 

Given  a set  of  such  measurements  {xj}™*  and  associated  index  functions  ( j(i),k(i )) 
specifying  the  targets  and  artefact  positions,  estimates  of  the  model  parameters  can  be 
determined  by  solving  a non-linear  least  squares  problem 


min 

{yjbftfcl.b 


XXfi> 


i=l 


where  f^y,-^),  b)  = x;  + e(xj,  b)  - yjtk. 

The  model  involves  three  sets  of  the  parameters:  the  target  locations  {y^},  transform- 
ation parameters  {t*,}  and  the  error  parameters  b.  Each  observation  equation  depends 
on  only  one  target  and  one  transformation,  so  that  the  Jacobian  matrix  J of  partial 
derivatives  can  be  ordered  to  have  a block-angular  structure  [2] 

' K 1 Ja  ‘ 

K2  J2 

J=  : « 

Kmx  Jmx  _ 

where  Kj  corresponds  to  the  parameters  y j and  the  border  blocks  {Jj}  correspond  to 
the  border  parameters  a = {{t*,},  b}.  The  frame  of  reference  for  the  targets  {y.,}  can  be 
specified  by  applying  three  appropriate  linear  equality  constraints  on  the  transformation 
parameters  {tfc}. 

While  scale  and  orthogonality  errors  are  often  major  contributors  to  the  systematic 
error  behaviour  of  a CMM,  there  is  no  guarantee  nor  does’ experience  show  that  they 
explain  the  full  extent  of  the  behaviour.  For  this  reason,  more  comprehensive  models  have 
been  developed  [3,  9].  However,  they  all  depend  on  the  approximation  of  actual  behaviour 
by  empirical  functions  such  as  polynomials  and  the  adequacy  of  the  approximation  is 
often  difficult  and  expensive  to  evaluate.  However,  if  we  always  rotate  and  translate  the 
artefact  according  to  the  symmetries  of  the  reference  artefact  so  that  the  targets  are 
always  located  (nominally)  at  a subset  of  a fixed  grid  of  points  in  the  CMM’s  working 
volume,  then  measurements  are  made  at  a finite  number  of  machine  locations.  To  the 
Zth  location  we  associate  a machine  error  e;.  If  the  ith  measurement  is  made  at  the  Ztli 
location  then  the  observation  equation  corresponding  to  (3.1)  is 

Xj  + e*  = yjjfc  + Cj. 

The  advantage  of  this  error  model  is  that  it  entails  no  significant  approximation:  the 
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Fig.  1.  Sparsity  structure  of  the  transpose  of  the  Jacobian  matrix  associated  with  the 
measurement  of  a 5 x 5 hole  plate  in  eight  positions. 

systematic  errors  are  modelled  exactly.  An  apparent  disadvantage  is  that  there  are  likely 
to  be  as  many  error  parameters  as  target  parameters  giving  rise  to  a sparsity  structure 
in  the  Jacobian  matrix  for  which  direct,  structure-exploiting  methods  provide  relatively 
minor  efficiency  gains.  Figure  1 shows  on  the  left  the  sparsity  structure  of  the  Jacobian 
matrix  J associated  with  the  measurement  of  a 5 x 5 hole  plate  in  eight  positions,  the 
first  four  corresponding  to  rotations  by  0,  90,  180  and  270  degrees,  the  second  four 
incorporating  a translation  as  well  as  a rotation.  In  each  position  the  location  of  the 
targets  y7-  are  measured  in  order.  The  nonzero  elements  of  the  matrix  are  represented 
by  a dot.  The  first  (second)  50  columns  correspond  to  the  derivatives  with  respect  to 
the  machine  error  parameters  e;  (target  parameters  y^)  and  the  last  24  correspond  to 
the  eight  sets  of  transformation  parameters  t^.  On  the  right  the  sparsity  structure  of  the 
triangular  factor  of  J is  illustrated  and  shows  the  substantial  fill-in  that  occurs. 

In  the  next  two  sections,  we  describe  approaches  for  dealing  efficiently  with  block- 
angular  and  more  general  sparse-block  structure. 

4 Algorithms  for  block-angular  systems 

We  consider  non-linear  least  squares  problems  where  the  optimisation  parameters  can 
be  partitioned  into  two  sets  r]  = fy,}"1'  and  a,  and  such  that  each  observation  equation 
involves  a and  at  most  one  set  of  parameters  yj.  Corresponding  to  (2.1),  we  have  instead 
an  objective  function  of  the  form 

F(?7,a)  = f0T(a)f0(a)  + 5^f7(yj!a)f,(yj,a). 

j 
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The  associated  Jacobian  matrix  J and  its  triangular  factor  R can  be  arranged  to  have 
the  form 


' Kx 

Ji  ' 

" Ri 

Bx 

J = 

k2 

Ji 

, R = 

R-2 

b2 

K-ny  Jny 
Jo 

J^ny  Buy 

Bo  , 

The  nonzero  blocks  of  the  matrix  R can  be  stored  compactly  in  a vector  r,  row  by  row. 

Efficient  updating  strategies  for  such  triangular  factors  have  been  incorporated  into 
a non-linear  least-squares  solver  to  deal  with  block-angular  problems.  It  is  assumed  that 
the  Jacobian  matrix  is  composed  of  ns  blocks  of  rows,  with  the  ith  block  depending  on 
at  most  one  set  of  parameters  yj,  j = j(i).  The  user  is  required  to  supply  a function  and 
gradient  evaluation  module  that  given  tj,  a and  1 < i < ns,  returns  j = j(i)  and 

f»(a),  Ju  3-  0, 

fi(yj.a),  Ju  Ku  j > 0. 

For  each  i,  the  triangular  factor  and  righthand  side  vector  is  updated  by  the  ith  block 
of  rows:  _ „ _ 


Rj(i) 

Bj(i) 

1 — > 

1 

I3 

e*. 

1 

Rq 

1 — > 

Ro 

Ri 

Ji 

L 0 Ji  j 

. Ji  . 

0 

Linear  equality  constraints  on  the  border  parameters  a implemented  using  the  orthogonal 
projection  approach  can  be  incorporated  by  setting  A :=  J»V 2 at  the  appropriate  stage. 

5 Algorithms  for  sparse-block  matrices 

Let  m x n matrix  S be  composed  of  ns  submatrices  Sk  of  dimension  mk  x nk-  We 
assume  that  Sk  is  stored  (column- wise  or  row- wise)  as  a column  vector  s*,.  The  inform- 
ation in  S can  be  encoded  in  a column  vector  s j and  an  indexing  set  Is  such  that 
7s(l  : 5 ,k)  = {ik-jk,rnk,nk-,  Ik)  where  ( ik,jk ) specifies  the  location  of  5^(1, 1)  in  S and 
Ik  indicates  that  sk  = Si(lk  ■ h + mknk  - 1)-  Blocks  of  such  matrices  can  be  easily 
represented  by  concatenating  the  s-vectors  and  index  matrices  Is  and  performing  some 
trivial  index  modifications.  Matrix- vector  multiplications  of  the  form  y :=  aSx  + fiy  are 
easily  implemented  through  a sequence  of  full  matrix  multiplications:  y :=  /?y,  followed 
by 

y(ik  : ik  + mk  - 1)  :=  y (ik  : h + mk  - 1)  + a5fcx(jfc  : jk  + nk  - 1), 

k = 1, . . . ,riB-  A similar  scheme  calculates  x :=  aST y + /?x.  The  storage  and  multiplica- 
tion scheme  can  be  modified  to  take  into  account  the  type  or  structure  of  the  submatrices 

Sk- 

To  implement  linear  equality  constraints,  it  is  required  to  perform  matrix  multiplic- 
ation by  a submatrix  V2  of  the  orthogonal  factor  of  the  constraint  matrix  C.  A simple 
scheme  can  be  implemented  using  the  LAPACK  routines  DGEQRF  (orthogonal  fac- 
torisation) and  DORMQR  (matrix  multiplication  by  an  orthogonal  matrix  stored  as  a 
product  of  Householder  matrices)  [8] . 
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FlG.  2.  Residual  errors  associated  with  the  first  1000  observations  for  models  a)  with 
no  error  separation  (dots)  and  b)  with  error  separation. 

We  have  implemented  a non-linear  least  squares  solver  for  sparse-block  systems.  The 
user  is  required  to  supply  a module  that  takes  as  input  the  current  estimate  a of  the 
optimisation  parameters  and  outputs  the  function  values  f(a)  and  the  Jacobian  matrix 
stored  in  sparse-block  form  (sj,Is).  The  solver  implements  a Gauss-Newton  approach 
using  the  LSQR  solver  to  find  the  Gauss-Newton  step  and  caters  in  a straightforward 
way  for  linear  equality  constraints.  The  solver  has  been  successfully  tested  in  a number 
of  self-calibration  problems.  For  example,  it  was  used  recently  in  the  calibration  of  a 
13  x 13  grid  of  targets  on  a glass  plate  by  a CMM  with  an  optical  probing  system.  The 
problem  involved  approximately  15,000  observation  equations  in  over  800  optimisation 
parameters  and  was  solved  in  a few  tens  of  seconds  using  a standard  laboratory  PC  (450 
MHz).  The  advantage  of  the  error  separation  model  is  illustrated  in  Figure  2 which  shows 
the  residual  errors  associated  with  the  first  1000  observations  for  models  a)  with  no  error 
separation  (dots)  and  b)  with  error  separation.  The  fit  for  the  error  separation  model  is 
much  superior.  The  practical  metrological  consequence  of  adopting  the  enhanced  model 
is  that  uncertainties  associated  with  the  target  locations  can  be  reduced  by  a factor 
of  five.  Importantly,  because  the  model  is  a realistic  approximation  of  the  measuring 
system,  we  can  have  confidence  in  the  uncertainty  estimates  derived  from  the  model. 

6 Concluding  remarks 

The  move  to  more  accurate  measurement  systems  has  led  to  more  comprehensive  models 
of  the  measuring  instrument  and  its  interaction  with  the  physical  quantity  being  meas- 
ured. These  models  include  parameters  that  describe  properties  of  the  instrument  and 
those  of  the  measurand.  The  aim  of  self-calibration  experiments  is  to  determine  as  much 
as  possible  about  both  sets  of  parameters  from  a set  of  measurement  experiments.  For 
models  with  a small  to  modest  set  of  parameters,  a full  matrix  approach  may  be  accept- 
able. For  larger  systems,  exploitation  of  sparsity  structure  in  the  defining  equations  is 
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highly  desirable  and  often  a stark  necessity  if  the  computations  are  to  be  made  in  an  ac- 
ceptable time  using  the  computing  resources  to  hand.  The  exploitation  of  block-angular 
structure  has  been  well-known  and  well-used  in  some  areas  of  metrology.  The  supporting 
numerical  technology  based  on  structured  orthogonal  factorisations  is  mature,  compact 
and  easily  implemented  using  standard  numerical  linear  algebra.  However,  these  tech- 
niques could  be  applied  more  widely  in  metrology,  making  feasible  approaches  that  have 
to  be  rejected  if  full  matrix  methods  only  are  to  be  used. 

The  use  of  sparse  matrix  techniques  is  relatively  rare  within  metrology.  We  have 
attempted  to  show  here  that  in  self-calibration  problems  in  dimensional  metrology,  they 
allow  us  to  develop  improved  models  that  provide  vastly  superior  fits  to  the  data,  with 
corresponding  improvements  in  the  evaluated  uncertainties  in  the  fitted  parameters.  The 
supporting  numerical  technology  is  maturing  and  accessible. 
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